####Annotations####
# File Name:    replication_jpp_limberg.R
# Author:       Julian Limberg
# Purpose:      Replication material for "What's Fair? Preferences for Tax Progressivity in the Wake of the Financial Crisis" (Journal of Public Policy)
#Last run:      R version 3.5.0 (2018-04-23) on 2 December 2018

####1. Load Packages####

#Install packages from CRAN
p_needed <- c("tidyverse","magrittr", "plyr", "ggplot2", "dplyr", "lme4", "ordinal", "texreg", "sandwich", "car", "grid", "stringr", "interplot","Hmisc", "xtable", "ggvis", "stargazer","simpleboot","margins","devtools","rlme","ggrepel","estimatr","lmtest","gridExtra","cowplot","usdm","extrafont","usdm")
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}

lapply(p_needed, require, character.only = TRUE)

#Set Seed
set.seed(2789)

####2. Load Data####

#Set working directory to current folder
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
data_2009 <- readRDS(file = "data_2009.rds")
data_1999_2009 <- readRDS(file = "data_1999_2009.rds")
datapit_2000_2016 <- readRDS(file = "pit_2000-2016.rds")
growth_1990_2014 <- readRDS(file = "growth_1990_2014.rds")

#Calculate aggregate change 1999 - 2009 (country weights)
data_1999_2009$taxprog_weight <- data_1999_2009$taxprog*data_1999_2009$WEIGHT
data_1999_2009$age_weight <- data_1999_2009$age*data_1999_2009$WEIGHT
data_1999_2009$parttime_weight <- data_1999_2009$parttime*data_1999_2009$WEIGHT
data_1999_2009$religious_weight <- data_1999_2009$religious*data_1999_2009$WEIGHT
data_1999_2009$unempl_weight <- data_1999_2009$unempl*data_1999_2009$WEIGHT

mean_1999_2009 <- 
  data_1999_2009 %>%
  dplyr::group_by(country,year) %>%
  dplyr::summarise(taxprog_weight = mean(taxprog_weight, na.rm = T),
                   age_weight = mean(age_weight, na.rm = T),
                   parttime_weight = mean(parttime_weight, na.rm = T),
                   religious_weight = mean(religious_weight, na.rm = T),
                   unempl_weight = mean(unempl_weight, na.rm = T),
                   wdi_gdpgr = mean(wdi_gdpgr, na.rm = T))

diff_1999_2009  <- mean_1999_2009 %>%
  dplyr::group_by(country) %>%
  dplyr::mutate(d_taxprog_weight = taxprog_weight - dplyr::lag(taxprog_weight), 
                d_age_weight = age_weight - dplyr::lag(age_weight), 
                d_parttime_weight = parttime_weight - dplyr::lag(parttime_weight), 
                d_religious_weight = religious_weight - dplyr::lag(religious_weight), 
                d_unempl_weight = unempl_weight - dplyr::lag(unempl_weight),
                dummy_growth = case_when(wdi_gdpgr < -2 ~ 1,
                                         TRUE ~ 0),
                dummy_growth1 = case_when(wdi_gdpgr < -1 ~ 1,
                                          TRUE ~ 0),
                dummy_growth2 = case_when(wdi_gdpgr < 0 ~ 1,
                                          TRUE ~ 0)) %>%
  dplyr::filter(year == 2009)

diff_1999_2009 <- data.frame(diff_1999_2009)
diff_1999_2009 <- separate(diff_1999_2009, country, c("country_short", "country_full"), sep = "-")
diff_1999_2009$country_short[diff_1999_2009$country_short == "United Kingdom (Great Britain only)"] <- "UK"
diff_1999_2009$country_full[diff_1999_2009$country_short == "UK"] <- "UK"
diff_1999_2009$country_full[diff_1999_2009$country_full == "Slovak Republic"] <- "Slovakia"
diff_1999_2009$country_full <- as.character(diff_1999_2009$country_full)

row.names(diff_1999_2009) <- diff_1999_2009$country_full

####3. Summary Statistics####

#Table OA2
summary_stats <- data_2009 %>%
  dplyr::group_by(country) %>%
  dplyr::summarise(taxprog_mean = mean(taxprog),
                   taxprog_sd = sd(taxprog),
                   importance_family_mean = mean(importance_family,na.rm=T),
                   importance_family_sd = sd(importance_family, na.rm=T),
                   payment_good_performance_mean = mean(payment_good_performance,na.rm=T),
                   payment_good_performance_sd = sd(payment_good_performance, na.rm=T),
                   corrupt_mean = mean(corrupt,na.rm=T),
                   corrupt_sd = sd(corrupt, na.rm=T),
                   growth2009 = mean(growth2009)) %>%
  arrange(desc(taxprog_mean))
summary_stats <- setNames(summary_stats, c("Country","Tax Prog. (Mean)","Tax Prog. (SD)","Des. Backgr. (Mean)","Des. Backgr. (SD)","Des. Behav. (Mean)","Des. Behav. (SD)","Des. Inst. (Mean)","Des. Inst. (SD)","Growth 2009"))
print(xtable(summary_stats,
             caption="Summary Statistics of Dependant Variables and the Main Macro Variable",
             digits=2))

#Table OA3
#Micro
micro_vars <- subset(data_2009, select = c(taxprog, importance_family, payment_good_performance, corrupt, parttime, unempl, educ, retired, nonlf, degree, age, sex, religious, dec_income))
stargazer(data.frame(micro_vars), type = "latex", summary = TRUE, covariate.labels=c("Tax Progressivity","Des. Backgr.","Des. Behav.","Des. Inst.","Part-Time Employed","Unemployed","In Education","Retired","Not in Labour Force","Educational Level","Age","Male","Religiosity","Income Deciles"))

#Macro
macro_vars <- data_2009 %>%
  dplyr::group_by(country) %>%
  dplyr::summarise(growth2009 = mean(growth2009),
                   growth2007 = mean(growth2007),
                   zscore = mean(zscore),
                   ln_gdp2009 = mean(ln_gdp2009),
                   gini_net = mean(gini_net),
                   gini_market = mean(gini_market),
                   introduction_pit = mean(introduction_pit),
                   sales_rev = mean(sales_rev),
                   soc_ben_2009 = mean(soc_ben_2009))

stargazer(data.frame(macro_vars), type = "latex", summary = TRUE, covariate.labels=c("Growth 2009","Growth 2007","Z-Score","GDP 2009 (ln)","Net Gini","Market Gini","Introduction PIT","Revenue from Sales Tax","Social Benefit Expenditure"))


####4. Summary Graphs####
#Figure 1
figure_1_growth <- growth_1990_2014 %>%
  ggplot(., aes(x=year, y=growth)) +
  geom_line() +
  labs(x = "Year", 
       y = "Real GDP Growth Rate") +
  theme_bw(base_size = 17) +
  theme_bw() +
  scale_x_continuous(limits = c(1990,2014), expand = c(0.01,0.01)) +
  theme(text = element_text(size=25, family = "Sabon"),
        panel.grid = element_blank())

#ggsave("../03_graphs/Fig_1_growth.pdf", plot = figure_1_growth, width = 22, height = 22, units = "cm", dpi = 800)

#Figure 3
pit_graph <- datapit_2000_2016 %>%
  dplyr::group_by(year) %>%
  dplyr::summarise(top_pit = mean(top_pit)) %>%
  dplyr::mutate(before = case_when(year < 2008 ~ 0,
                                   TRUE ~ 1))

figure_3_pit_rates <- pit_graph %>%
  ggplot(., aes(x=year, y=top_pit)) +
  geom_line(size = 1) +
  xlab("Year") +
  ylab("Top PIT Rate") +
  coord_cartesian(xlim = c(2000,2016), ylim = c(40,50)) +
  theme_bw() +
  theme(text = element_text(size=25, family = "Sabon"),
        panel.grid = element_blank())

#ggsave("../03_graphs/Fig_3_rates.pdf", plot = figure_3_pit_rates, width = 22, height = 22, units = "cm", dpi = 800)

#Figure OA1
summary_graph <- data_2009 %>%
  dplyr::group_by(country) %>%
  dplyr::mutate(taxprog_weight = taxprog * WEIGHT) %>%
  dplyr::summarise(taxprog_mean = mean(taxprog_weight),
                   growth2009 = mean(growth2009))

#Figure OA1 (a) - full sample
figure_A1_full <- summary_graph %>%
  ggplot(.,aes(x=growth2009,y=taxprog_mean)) +
  geom_point() +
  ggrepel::geom_text_repel(label=summary_graph$country, 
                           data=summary_graph,
                           force=3, 
                           size = 6,
                           min.segment.length=2, 
                           family = "Sabon") +
  geom_smooth(method = lm, colour = "black", se = F) +
  xlab("Growth 2009") +
  ylab("Preferences for Tax Progressivity") +
  theme_bw() +
  theme(text = element_text(size=25, family = "Sabon"),
        panel.grid = element_blank())

#ggsave("../03_graphs/Fig_OA1_scatterplot_a.pdf", plot = figure_A1_full, width = 22, height = 22, units = "cm", dpi = 800)

#Figure OA1 (b) - excluding outliers
summary_graph_reduced <- summary_graph %>%
  dplyr::filter(!(country %in% c("Ukraine","Latvia","Estonia","South Korea","Russia"))) 

figure_A1_subset <- summary_graph_reduced %>%
  ggplot(.,aes(x=growth2009,y=taxprog_mean)) +
  geom_point() +
  ggrepel::geom_text_repel(label=summary_graph_reduced$country, 
                           data=summary_graph_reduced,
                           force=3, 
                           size = 6,
                           min.segment.length=2, 
                           family = "Sabon") +
  geom_smooth(method = lm, colour = "black", se = F) +
  xlab("Growth 2009") +
  ylab("Preferences for Tax Progressivity") +
  theme_bw() +
  theme(text = element_text(size=25, family = "Sabon"),
        panel.grid = element_blank())

#ggsave("../03_graphs/Fig_OA1_scatterplot_b.pdf", plot = figure_A1_subset, width = 22, height = 22, units = "cm", dpi = 800)

####5. Main Analyses####

#Table 1
#1: Bivariate multilevelmodel without micro variables
summary(mod1 <- lmer(taxprog ~ growth2009 + (1|country), data = data_2009))

#2: Multilevel models with micro controls
summary(mod2 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + growth2009 + (dec_income|country), data = data_2009, REML = FALSE))

#3: Growth for first half of 2009 for countries where survey has taken place earlier
summary(mod3 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income+ importance_family + corrupt + payment_good_performance + growth091 + (dec_income|country), data = data_2009, REML = FALSE))

#4: Add macro controls
summary(mod4 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + growth2009 + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

screenreg(list(mod1,mod2,mod3,mod4), stars = c(0.01, 0.05, 0.10) ,  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, dcolumn = TRUE, longtable = FALSE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models for Tax Progressivity", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.","growth2009" = "Growth 2009", "growth091" = "Growth First Half 2009", "growth2007" = "Growth 2007", "zscore" = "Z-Score", "ln_gdp2009" = "GDP 2009 (ln)","(Intercept)" = "(Intercept)"))

#Multicollinearity of the Micro variables
df <- data_2009 %>%
  dplyr::select(parttime,unempl,educ,retired,nonlf,degree,age,sex,religious,dec_income,importance_family,payment_good_performance,corrupt) %>%
  data.frame()
vif(df)

res <- cor(df)
round(res, 2) %>% View()

#Table 2
#1: Add pre-crisis government size
summary(mod1.0 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + zscore + ln_gdp2009 + growth2007 + imf_exp + (dec_income|country), data = data_2009, REML = FALSE))

#2: Output gap 2009
summary(mod1.1 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income+ importance_family + corrupt + payment_good_performance + output_gap + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

#3: Output gap 2009/2010
summary(mod1.2 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + output_gap_2 + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

#4: Check with 1999 - placebo
dat_1999 <- data_1999_2009 %>% 
  dplyr::filter(year==1999) 
summary(mod1.3 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + payment_good_performance + corrupt + wdi_gdpgr + growth1997 + zscore + ln_wdi_gdpcapcur99 + (dec_income|country), data = dat_1999, REML = FALSE))

#5: Comparison 1999
summary(mod1.4 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + payment_good_performance + corrupt + growth1999 + growth1997 + zscore + ln_wdi_gdpcapcur99 + (dec_income|country), data = dat_1999, REML = FALSE))

#6: Same analyses for 2009 with reduced country sample
data_2009_red <- data_2009 %>%
  dplyr::filter(country %in% c("Australia","United Kingdom","Austria","Chile","France","Germany","Hungary","Japan","Latvia","New Zealand","Norway","Philippines","Poland","Russia","Slovak Republic","Slovenia","Spain","Sweden","United States")) %>%
  dplyr::rename(sexMale = sex,
                wdi_gdpgr = growth2009)

summary(mod1.5 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sexMale + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + wdi_gdpgr + growth2007 + zscore + ln_gdp2009 + (dec_income|country), data = data_2009_red, REML = FALSE))

screenreg(list(mod1.0,mod1.1,mod1.2,mod1.3,mod1.4,mod1.5), stars = c(0.01, 0.05, 0.10) ,  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, dcolumn = TRUE, longtable = FALSE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models for Tax Progressivity in 2009 and 1999", caption.above = TRUE, custom.model.names = c("Control Exp.","Output Gap I","Output Gap II", "Placebo 1999", "Reduced 1999", "Reduced 2009") , custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "sexMale" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.", "growth2009" = "Growth 2009", "wdi_gdpgr" = "Growth 2009", "growth1999" = "Growth 1999", "output_gap" = "Output Gap 2009",  "output_gap_2" = "Output Gap 2009/2010", "growth2007" = "Growth t-2", "growth1997" = "Growth t-2", "zscore" = "Z-Score", "ln_gdp2009" = "GDP (ln)", "ln_wdi_gdpcapcur99" = "GDP (ln)", "imf_exp" = "Government Exp.","(Intercept)" = "(Intercept)"))

#Table 3
#1
summary(mod_diff1 <- lm_robust(d_taxprog_weight ~ wdi_gdpgr, se_type = "HC1", data = diff_1999_2009))

#2
summary(mod_diff2 <- lm_robust(d_taxprog_weight ~ wdi_gdpgr + d_age_weight, se_type = "HC1", data = diff_1999_2009))

#3
summary(mod_diff3 <- lm_robust(d_taxprog_weight ~ wdi_gdpgr + d_age_weight + d_parttime_weight, se_type = "HC1", data = diff_1999_2009))

#4
summary(mod_diff4 <- lm_robust(d_taxprog_weight ~ wdi_gdpgr + d_age_weight + d_parttime_weight + d_religious_weight, se_type = "HC1", data = diff_1999_2009))

#5
summary(mod_diff5 <- lm_robust(d_taxprog_weight ~ wdi_gdpgr + d_age_weight + d_parttime_weight + d_religious_weight + d_parttime_weight + d_unempl_weight, se_type = "HC1", data = diff_1999_2009))

#6
summary(mod_diff6 <- lm_robust(d_taxprog_weight ~ dummy_growth + d_age_weight + d_parttime_weight + d_religious_weight + d_parttime_weight + d_unempl_weight, se_type = "HC1", data = diff_1999_2009))

screenreg(list(mod_diff1,mod_diff2,mod_diff3,mod_diff4,mod_diff5,mod_diff6), 
       stars = c(0.01, 0.05, 0.10), 
       no.margin = TRUE, 
       include.adjrs = FALSE,
       include.rmse = FALSE,
       booktabs = TRUE, 
       dcolumn = TRUE, 
       longtable = FALSE, 
       center = TRUE, 
       digits = 4, 
       caption = "Determinants of Change in Preferences for Tax Progressivity 1999--2009",
       caption.above = TRUE, 
       custom.coef.map = list("wdi_gdpgr" = "Growth 2009", "dummy_growth" = "Major Economic Crisis", "d_age_weight" = "$\\Delta$ Age", "d_parttime_weight" = "$\\Delta$ Part-Time Empl.", "d_religious_weight" = "$\\Delta$ Religosity", "d_unempl_weight" = "$\\Delta$ Unempl.","(Intercept)" = "(Intercept)"))


####6. Cross-Level Interaction Effects####

#Table OA5
summary(mod_clie_1 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family * growth2009 + corrupt + payment_good_performance + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

summary(mod_clie_2 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance * growth2009 + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

summary(mod_clie_3 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt * growth2009 + payment_good_performance + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_2009, REML = FALSE))

screenreg(list(mod_clie_1,mod_clie_2,mod_clie_3), stars = c(0.01, 0.05, 0.10) ,  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models with Interaction Effects", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.","growth2009" = "Growth 2009", "importance_family:growth2009" = "Des. Backgr * Growth 2009", "payment_good_performance:growth2009" = "Des. Behav. * Growth 2009", "corrupt:growth2009" = "Des. Inst. * Growth 2009", "growth2007" = "Growth 2007", "zscore" = "Z-Score", "ln_gdp2009" = "GDP 2009 (ln)","(Intercept)" = "(Intercept)"))

#Number countries
growth_country <- data_2009 %>%
  dplyr::group_by(country) %>%
  dplyr::summarise(growth = mean(growth2009))

hist_countries <- growth_country %>%
  ggplot(., aes(x=growth)) + 
  geom_histogram(col="black", 
                 alpha = 0) +
  coord_cartesian(xlim = c(-15,3)) +
  labs(x = "Growth 2009", 
       y = "No. of Countries") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        text = element_text(size = 25, family = "Sabon"))

#Interaction Graphs
p1 <- interplot(mod_clie_1, var1 = "importance_family", var2 = "growth2009", hist = FALSE, ci = 0.90) +
  coord_cartesian(xlim = c(-15,3),ylim = c(-0.02,0.08)) +
  labs(y = "Marg. Effect Des. Backgr.") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        text = element_text(size = 25, family = "Sabon")) +
  geom_hline(aes(yintercept = 0), linetype = "dashed")

figure_2_plot_1 <- plot_grid(p1, hist_countries, align = "v", nrow = 2,rel_heights = c(3/4, 1/4))

#ggsave("../03_graphs/Fig_2_interaction_1.pdf", plot = figure_2_plot_1, width = 22, height = 22, units = "cm", dpi = 800)

p2 <- interplot(mod_clie_2, var1 = "payment_good_performance", var2 = "growth2009", hist = FALSE, ci = 0.90, var2_dt = summary_graph$growth2009) +
  coord_cartesian(xlim = c(-15,3),ylim = c(-0.02,0.08)) +
  labs(y = "Marg. Effect Des. Behav.")+
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        text = element_text(size = 25, family = "Sabon")) +
  geom_hline(aes(yintercept = 0), linetype = "dashed")

figure_2_plot_2 <- plot_grid(p2, hist_countries, align = "v", nrow = 2,rel_heights = c(3/4, 1/4))

#ggsave("../03_graphs/Fig_2_interaction_2.pdf", plot = figure_2_plot_2, width = 22, height = 22, units = "cm", dpi = 800)

p3 <- interplot(mod_clie_3, var1 = "corrupt", var2 = "growth2009", hist = FALSE, ci = 0.90, var2_dt = summary_graph$growth2009) +
  coord_cartesian(xlim = c(-15,3),ylim = c(-0.02,0.08)) +
  labs(y = "Marg. Effect Des. Inst.") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        text = element_text(size = 25, family = "Sabon")) +
  geom_hline(aes(yintercept = 0), linetype = "dashed")

figure_2_plot_3 <- plot_grid(p3, hist_countries, align = "v", nrow = 2,rel_heights = c(3/4, 1/4))

#ggsave("../03_graphs/Fig_2_interaction_3.pdf", plot = figure_2_plot_3, width = 22, height = 22, units = "cm", dpi = 800)

####7. Additional Robustness Checks####

#Table OA4: Additional Macro IVs
summary(mod2.6 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + gini_net + (dec_income|country), data = data_2009, REML = FALSE))
summary(mod2.7 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + gini_market + (dec_income|country), data = data_2009, REML = FALSE))
summary(mod2.8 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + introduction_pit + (dec_income|country), data = data_2009, REML = FALSE))
summary(mod2.9 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + sales_rev + (dec_income|country), data = data_2009, REML = FALSE))
summary(mod2.10 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + left + center + no_affil + importance_family + corrupt + payment_good_performance + growth2009 + soc_ben_2009 + (dec_income|country), data = data_2009, REML = FALSE))

screenreg(list(mod2.6,mod2.7,mod2.8,mod2.9,mod2.10), stars = c(0.01, 0.05, 0.10),  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 3, no.margin =  TRUE, caption = "Results Multilevel Models for Tax Progressivity: Robustness Checks", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.", "corrupt" = "Des. Inst.", "growth2009" = "Growth 2009","gini_net" = "Net Gini", "gini_market" = "Market Gini", "introduction_pit" = "Introduction PIT", "sales_rev" = "Revenue from Sales Taxes", "soc_ben_2009" = "Social Benefits Expenditure","(Intercept)" = "(Intercept)"))

#Table OA6: GLM
summary(mod2.0o <- clmm(factor(taxprog) ~ growth2009 + (1|country), data = data_2009))
summary(mod2.1o <- clmm(factor(taxprog) ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + growth2009 + (1|country), data = data_2009))
summary(mod2.2o <- clmm(factor(taxprog) ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income+ importance_family + corrupt + payment_good_performance + growth091 + (1|country), data = data_2009))
summary(mod2.3o <- clmm(factor(taxprog) ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + growth2009 + zscore + ln_gdp2009 + growth2007 + (1|country), data = data_2009))

screenreg(list(mod2.0o,mod2.1o,mod2.2o,mod2.3o), stars = c(0.01, 0.05, 0.10) ,  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models for Tax Progressivity", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.","growth2009" = "Growth 2009", "growth091" = "Growth First Half 2009", "growth2007" = "Growth 2007", "zscore" = "Z-Score", "ln_gdp2009" = "GDP 2009 (ln)","(Intercept)" = "(Intercept)"))

#Table OA7: Additional Political Affiliation Variables
data_political <- data_2009 %>% 
  filter(!(country %in% c("Chile", "Hungary", "Israel"))) #No data on political affiliation
summary(mod2.12 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + left + center + no_affil + growth2009 + (dec_income|country), data = data_political, REML = FALSE))

summary(mod2.13 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income+ importance_family + corrupt + payment_good_performance + left + center + no_affil + growth091 + (dec_income|country), data = data_political, REML = FALSE))

summary(mod2.14 <- lmer(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + left + center + no_affil + growth2009 + zscore + ln_gdp2009 + growth2007 + (dec_income|country), data = data_political, REML = FALSE))

screenreg(list(mod2.12,mod2.13,mod2.14), stars = c(0.01, 0.05, 0.10),  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, include.ci = FALSE, dcolumn = TRUE, longtable = TRUE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models for Tax Progressivity: Additional Party Variables", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.", "left" = "Left Affiliation","center" = "Center Affiliation" , "no_affil" = "No Affiliation","growth2009" = "Growth 2009", "growth091" = "Growth First Half 2009", "growth2007" = "Growth 2007", "zscore" = "Z-Score", "ln_gdp2009" = "GDP 2009 (ln)","(Intercept)" = "(Intercept)"))

#Table OA8: Robustness Micro Variables
summary(mod2.10 <- lm(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance + factor(country), data = data_2009)) #Country FE

summary(mod2.11 <- lm_robust(taxprog ~ parttime + unempl + educ + retired + nonlf + degree + age + sex + religious + dec_income + importance_family + corrupt + payment_good_performance, clusters = country, se_type = "stata", data = data_2009)) #Country clustered SE

screenreg(list(mod2.10,mod2.11), stars = c(0.01, 0.05, 0.10) ,  include.bic = FALSE, include.variance = FALSE, booktabs = TRUE, include.ci = FALSE, dcolumn = TRUE, longtable = TRUE, center = TRUE, no.margin =  TRUE, digits = 4, caption = "Results Multilevel Models for Tax Progressivity: Robustness Checks", caption.above = TRUE, custom.coef.map = list("parttime" = "Part-Time Employed", "unempl" = "Unemployed", "educ" = "In Education", "retired" = "Retired", "nonlf" = "Not in Labour Force", "degree" = "Educational Level", "age" = "Age", "sex" = "Male", "religious" = "Religiosity", "dec_income" = "Income", "importance_family" = "Des. Backgr.", "payment_good_performance" = "Des. Behav.",  "corrupt" = "Des. Inst.", "left" = "Left Affiliation","center" = "Center Affiliation" , "no_affil" = "No Affiliation","growth2009" = "Growth 2009", "growth091" = "Growth First Half 2009", "growth2007" = "Growth 2007", "zscore" = "Z-Score", "ln_gdp2009" = "GDP 2009 (ln)","(Intercept)" = "(Intercept)"))